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We present a method for computing parameter sensitivities and response coefficients in Brownian 
dynamics simulations. The method involves tracking auxiliary variables (Malliavin weights) in 
addition to the usual particle positions, in an unperturbed simulation. The Malliavin weights 
sample the derivatives of the probability density with respect to the parameters of interest and are 
also interesting dynamical objects in themselves. Malliavin weight sampling is simple to implement, 
applies to equilibrium or nonequilibrium, steady state or time-dependent systems, and scales more 
efficiently than standard finite difference methods. 
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The response of a system to infinitesimal changes in 
an external field provides important insights into its un- 
derlying physics. Divergences in such response func- 
tions can indicate the presence of phase transitions, while 
their relation to fluctuations in unperturbed systems via 
fluctuation-dissipation theorems (FDTs) provides a key 
diagnostic of the difference between equilibrium and non- 
equilibrium systems [l| . Knowledge of the response of a 
system to changes in internal parameters (e. g. those con- 
trolling inter-particle forces) is also of great importance, 
since it has the potential greatly to accelerate the fitting 
of force fields to experimental data. For equilibrium sys- 
tems, responses to perturbations can often be computed 
from the properties of unperturbed systems via known 
statistical mechanical relations However, for sys- 

tems which are far from steady-state, and/or whose dy- 
namics does not obey detailed balance, one is generally 
forced to resort to finite differencing: explicitly taking 
the difference between simulation trajectories generated 
at slightly different parameter values [5[ . While finite dif- 
ferencing can be made more efficient by reuse of random 
number streams [H, 0], one still has to re-simulate the 
perturbed system for each parameter of interest. 

In this Letter, we present a simple and generic method 
for computing responses to infinitesimal changes in in- 
ternal or external parameters in stochastic Brownian dy- 
namics simulations, which may be in or out of equilib- 
rium. The method does not require simulation of the 
perturbed system; instead, it involves tracking, in an un- 
perturbed system, auxiliary stochastic variables which 
sample the derivatives of the probability density with re- 
spect to the parameters of interest. We term these auxil- 
iary variables 'Malliavin weights' as the method has close 
links to the Malliavin calculus programme Q, used in 
quantitative finance for deriving price sensitivities (i. e. 
'Greeks') @. Our method extends approaches previously 
proposed for kinetic Monte-Carlo simulations [lol - [l^ | to 
a much wider set of problems. It also has interesting 



links to molecular dynamics methods in which response 
coefficients are computed via the integration of adjunct 
equations of motion for individual particles [H, [l(J • Since 
Brownian dynamics is very widely used (l3l | in the study 
of non-equilibrium statistical physics problems such as 
driven steady states, active soft matter, and modeling 
sub-cellular processes in biology, we anticipate that our 
method should prove widely applicable. 

We begin by considering a collection of i = 1 . . . N 
interacting particles undergoing overdamped Brownian 
motion, described by the coupled Langevin equations 



~dt 



Dfj 
k B T 



(1) 



Here r, and fi are the position of the ith particle and 
the force acting on it, respectively, D is the diffusion co- 
efficient (which for simplicity we here assume to be con- 
stant and the same for all particles), fee is Boltzmann's 
constant, T is temperature, and the ffi are independent 
vectors of Gaussian white noise of amplitude 2D. We 
now add an extra variable q\ (a Malliavin weight) which 
evolves according to 



dg\ _ 1 
dt ~ 2k B T 



v ^ 



(2) 



where A is a parameter of interest for which the only re- 
quirement is that dfi/dX is known. Note that q\ does 
not perturb the dynamics of the particles; it merely acts 
as a 'readout' and should be initialised to q\ = 0. The 
interpretation of Eq. ([2]) as a stochastic differential equa- 
tion is straightforward and uniquely defined — the practi- 
cal implementation is described in Supplementary Mate- 
rial 17). The noise vector ffi is identical in Eqs. ([1} and 
— in each Brownian dynamics timcstcp, q\ is updated 
using the same set of random numbers that were chosen 
in the update of the particle positions. Our central claim 
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is that for any function of the particle positions ^4({fi}) 



9(A) 
dX 



(3) 



i. e. the response of A to the parameter A is given by 
the average of A in the unperturbed system, weighted 
by the appropriate Malliavin weight q\. Eq. © has im- 
portant practical implications. Since the computation 
of q\ via Eq. ([2]) is independent of A, the same q\ can 
be used to compute the sensitivity of multiple system 
properties to the parameter A. Moreover, one can track 
multiple weights corresponding to different choices of A, 
with marginal additional cost. 

Eq. ([3]) is the key result of this Letter. It can be 
proved by taking moments of a Chapman-Kolmogorov 
equation for the evolution of the joint probability distri- 
bution P({r,}, q\; t) for the set of particle positions {fj} 
and the Malliavin weight q\. The details are given in 
Supplementary Material [l7|. A crucial intermediate 
result is that the conditional average of q\ for a given 
set of particle positions {ri}, which we denote (qx)^}, 
is given by 

/ \ _ Jdqxq\P({n},qx;t) <9 In P({f z }; 



JdqxP({n},q X ;t) 



dX 



where P({fi};t) is the probability distribution for the 
particle positions. Thus the Malliavin weight in fact sam- 
ples the conjugate 0] variable OhiP/dX. 

The proof makes no assumptions about the system be- 
ing in steady state or obeying detailed balance; thus our 
method is valid for systems far from steady state, or in 
driven steady states, as well as for those at equilibrium. 
Moreover, as we show in Supplementary Material [l7| . 
our approach can easily be extended to systems under- 
going underdamped Brownian motion and to the compu- 
tation of higher-order derivatives. 

To illustrate the method, we turn to a simple example 
for which analytical results are available: a single par- 
ticle in a one-dimensional harmonic trap described by a 
potential U = |ks 2 - hx. To make contact with linear re- 
sponse theory, we take the parameter of interest to be the 
strength of the applied external force h. We set the parti- 
cle mobility to unity so that the diffusion constant D = T 
where the temperature T is in units of fee- At equilib- 
rium, P cq {x) ~ e~ u / T , so that dlnP cq /dh = (x-(x))/T 
0, and the FDT holds: d(x)/dh = ((x 2 ) - (x) 2 )/T. 
We now apply Malliavin weight sampling to the time- 
dependent situation in which the particle starts from 
x = x at t — and relaxes towards its equilibrium 
position. Eqs. ((T|) and @ become 



dx 

lit 



-KX 



dqh 
dt 



2T 



(5) 




These equations can be solved exactly [17[ to give 



FIG. 1: (color online) Numerical simulation of Eqs. ((5J with 
ft = 2, x'o = 1, and D = T = 1. The falling curves show 
(x(i)) for h — (solid line) and h — 0.1 (dashed line). The 
rising curves show (x(t)qh(t)) (blue solid line) and d{x(t)) /dh 
evaluated using forward finite differencing (blue circles). Av- 
erages are over 10 replicate simulations. The inset shows the 
long-time behaviour of individual x(t) and qh(t) trajectories. 
While x(t) stays close to its equilibrium value (black), qh(t) 
(five replicates, red) executes a random walk. 



P(x,qh',t) as a bivariate Gaussian with 

(x) = x e~ Kt + (V«)(l ~ e " Kt ) , (Qh) = , 
(x 2 )-(x) 2 = (T/ K )(l-e- 2Kt ), (6) 
(xq h ) = (l-e-^)/ K , (q 2 )=t/2T. 

This result allows us to verify directly that Eq. (fj| is 
satisfied, that is to say (qh) x — Jdx qf l P(x,qh',t) — 
dhxP{x\t)/dh, where P(x;t) = Jdq h P(x,q h ;t). Fig. Q] 
shows simulation results for d(x) / dh, computed by Malli- 
avin weight sampling (MWS) and by forward finite dif- 
ferencing; the inset shows trajectories for q^ from repli- 
cate simulation runs. The Malliavin weight qh behaves 
as a random walk with zero mean and a diffusion co- 
efficient 1/4T. Eq. © shows that, by analogy with 
the equilibrium FDT, one can define a time-dependent 
effective temperature T cS /T = (1 - e~ 2Kt )/(l - e~ Kt ) 
such that d{x)/dh = ((x 2 ) — {x) 2 )/T e g. Calculating 
the conditional average of the Malliavin weight gives 
(qh)x = (x — (x))/T e g. Interestingly this has the same 
form as the equilibrium case but again features T G ff. 
Thus the effective temperature has a wider relevance than 
would be apparent from the time-dependent FDT since 
d(A)/dh = ((Ax) - (A)(x))/T eS , where A(x) is any func- 
tion of the particle position. 

An important practical issue is raised by the fact that 
qx(t) behaves as a random walk (inset to Figure QJ: to 
compute responses to parameter perturbations for sys- 
tems in steady state, we cannot simply monitor q\ for 
longer and longer times until the system reaches its 
steady state. This is because replicate trajectories of q\ 
diverge at long times and measurements of (A(t)q\(t)) 
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FIG. 2: (color online) Interacting particles in a two- 
dimensional harmonic trap under shear. The main plot shows 
how the correlation function (xyAq^) (solid line; calculated in 
an unsheared system) asymptotes to d{xy) /d^lj^o = 0.0890± 
0.0005 calculated by centered finite differencing (dashed line). 
The inset shows 100 superimposed simulation snapshots at 
At = 2, colored by individual particle Malliavin weight. 



incur a large sampling error. Fortunately, we can cir- 
cumvent this problem by computing instead the correla- 
tion function [13, El C(t,t ) = (A(t) Aq\(t, t )) where 
Aq\(t,to) = q\(t) — qx(to). In steady state this becomes 
a well-defined function of At = t — to and we expect it to 
obey C(At) -> d(A)/d\ as At -> oo [H[l3|. Like any 
other correlation function, C(t) converges to its asymp- 
totic value on a time scale set by the spectrum of relax- 
ation times in the problem. 

Next, we demonstrate the application of MWS to 
a much less trivial example: a non-equilibrium driven 
steady state formed by a cluster of particles in a two- 
dimensional harmonic trap, under shear. This example is 
motivated by recent experimental studies of colloidal par- 
ticles in optical traps, which have provided new insights 
into statistical physics at the microscale 15( . We suppose 
that the particles interact with each other via a repulsive 
screened Coulomb potential u(r) = (T /r)e~ r / Tc - , with 
coupling strength T and range r c . The total potential 



energy of the cluster is then U = J2i>j u ( r ij) + \ K T1 I 



(where k is the spring constant of the trap), and the 
Langevin equations for the particle positions (xi,yi) are 



dxi 
~dt 



dU 



+ iyi + Vx.i , 



dt 



dU_ 

dyi 



Vv,- 



(7) 



where 7 is the shear rate and r) x ,i and rj y .i are noise terms 
defined as in the previous example. We set r c = k^T = 
D = 1 to fix units of length, energy, and time, and 
choose N = 10, ft = 10 and T = 25; for this param- 
eter set, the particles form a dense, strongly correlated 
cluster in the trap (see inset to Fig. [21). To character- 
ize the morphology of the cluster we use quantities like 
(xy) = (J2i=i x iVi}/N- We first focus on the sensitivity 
of this quantity to changes in the shear rate 7. We there- 



fore track the Malliavin weight q^ which, from Eq. ([2]), 
obeys 



dq^j 1 
~dt = 2k B T^ yiVx '' 

i — 1 



(8) 



To compute d(xy)/d A f for the system in steady-state, 
we use the correlation function approach outlined 
above. Fig. [2] shows that, for 7 = 0, C(At) = 

<(l/JV)(Eiiafi(*)»i(*))(*r(*) - 9j(t - At))) tends to 
d(xy)/dj\j = o as At increases, where d(xy)/d , y\^-o (the 
dashed line) is calculated by finite differencing. In fact 
the Malliavin weight q^ turns out to be an interest- 
ing physical quantity in itself. By splitting the sum in 
Eq. (|5J) into individual particle contributions, one can 
track Malliavin weights q^ i for each individual particle. 
These provide insight into the response to shear of the 
one-particle probability distribution P(xi, j/j). As shown 
in the inset to Fig. [5p, the individual contributions to 
the Malliavin weight are biased towards being positive in 
the first and third quadrants and negative in the second 
and fourth quadrants, corresponding to the distortion of 
the particle cloud by the shear. 

An important feature of MWS is that, from a single 
simulation run, one can compute the sensivitities of any 
function of the particle coordinates, to any parameter 
of the system. One simply needs to track the Malliavin 
weights corresponding to all the parameters that are of 
interest, using the appropriate dynamical rules as derived 
from Eq. ([2]). For example in this problem q-p obeys 



dgr 
dt 



1 N 



d 2 U 



2k B T^\drdx< 

i—l 



d 2 U 

dTdy. 



■Vv, 



(9) 



and an analogous equation for q K is easily written down. 
Fig. [3^ shows the full panoply of responses of the cluster 
morphology parameters (xy) and (x 2 y 2 ) to the parame- 
ters of the problem, for T = 25, obtained from a single 
simulation run. Second order derivatives were computed 
as described in Supplementary Material [l7j ]. 

We now demonstrate how MWS can reveal subtle de- 
tails of how the response to shear, d(xy) /<97|^ = 0j depends 
on {x 2 y 2 } 1 ^ 2 , which provides a representative measure 
of the area of the particle cloud. For non-interacting 
particles one can obtain analytically [I?} the intriguing 
quasi-FDT result d{xy)/dj = {x 2 y 2 )/2T, which holds 
at r = 7 = and for all values of k. The case of in- 
teracting particles, however, cannot be solved analyt- 
ically. We therefore simulated the cloud at a series 
of increasing values of T (i. e. increasing cluster size) 
keeping k = 10. The results shown in Fig [3j) suggest 
that d(xy)/d A f\j = o is quite accurately proportional to 
(x 2 y 2 ) 1 / 2 . More generally we can define an effective ex- 
ponent (3 = d\n(d(xy)/d r y)/d\n.(x 2 y 2 ). This can be eval- 
uated as r varies at fixed k, or vice versa. In the former 
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FIG. 3: (color online) (a) Multiple parameter dependencies 
from a single simulation, (b) Dependence of shear response on 
the size of the cluster, as the interaction strength T is varied 
(the double circle is T = 0). The inset shows the effective 
exponent /3 computed from Eq. (|10|) in the main text, using 
MWS to evaluate the derivatives. 



case, expanding the derivatives gives 



d 2 (xy)/djdT 
d(x 2 y 2 )/dT d(xy)/dj 



(x 2 y 2 ) 



(10) 



We used MWS to compute the derivatives in Eq. (ITOl) . 
with results shown in the inset in Fig. [3]3. As expected 
from the main plot, f3 0.5, but a more subtle depen- 
dence is also apparent. An analogous calculation for the 
latter case (fixed T varying k) shows that the correspond- 
ing effective exponent increases from /3 = 1 at T = (the 
quasi-FDT result) to (3 « 1.3 at V = 25. Hence the inter- 
acting particle case shows considerably more complexity 
than the non-interacting one, and there is apparently lit- 
tle universality. 

These examples demonstrate that MWS provides a 
simple and easy-to-implcment alternative to finite dif- 
ferencing, for Brownian dynamics problems that may be 
time-dependent or in steady state, in or out of equilib- 
rium. Let us now discuss the question of efficiency. As 
highlighted above, in MWS one has access to response 
coefficients for all parameters of the problem (for which 
the derivatives dfi/dX in Eq. @ are known), with little 



additional cost, since integrating the equations of motion 
for the q\ requires no new random numbers and also typ- 
ically does not require recalculation of the forces. MWS 
also scales more efficiently with the computational effort 
than does standard finite differencing. For MWS, the 
error in a computation of d{A(t))/d\ scales as M~ 1//2 , 
where M is the number of replicate simulation runs used 
to compute the averages (this follows from the usual scal- 
ing of the standard error in the mean with the number 
of samples). For finite differencing, there is an inherent 
tradeoff between the systematic error introduced by us- 
ing a too-large perturbation and the random sampling 
error that arises when the perturbation is very small. 
One can show Q that the best possible choice of the per- 
turbation size results in an error that scales as i\/ _1 / 4 
for a forward finite differencing scheme, and M -1 / 3 for a 
centered scheme. Although the scaling can be improved 
to M' 1 / 2 by using a common random number scheme 
0, 0], even here we expect MWS to be more efficient, 
since it does not require the perturbed system to be ex- 
plicitly simulated. For example we note that more than 
six times as many force evaluations were used to calcu- 
late d(xy)/d A f in Fig. [2] by finite differencing, compared 
to calculating the same quantity to comparable accuracy 
by the MWS correlation function method. 

MWS has the potential greatly to facilitate the param- 
eterization of force fields, when combined with gradient- 
based search and optimisation algorithms. This should 
be especially relevant for mesoscale problems where 
Brownian dynamics algorithms such as dissipative parti- 
cle dynamics (DPD) are often the method of choice [3| 
(note that the application to DPD should take account of 
the fact that the noise terms are pairwise central random 
forces). An equally important application of MWS is in 
the computation of response functions to external fields — 
for example in the context of dynamical phase transi- 
tions. Interestingly, while we have focused here only on 
its long-time limit, the time-correlation function C(At) 
is actually equivalent to the time-dependent response to 
a step perturbation. This point will be explored in more 
detail in future work. A particularly interesting question 
concerns the extent to which the Malliavin weight itself 
can be used as an autonomous order parameter. For in- 
stance in the Id trap problem the conditional average 
(<lh)x features the time-dependent effective temperature 
T e s, generalising the linear response result. Certainly 
for glassy systems we expect that the correlation func- 
tions (AAq\) will show typical non-ergodic memory and 
aging phenomenology [10j: it is interesting to speculate 
whether the behaviour of q\ itself could be used as an 
alternative signature of glassy behaviour. 

We thank Mike Allen, Mike Cates, Alessandro Laio 
and Bartlomiej Waclaw for discussions. RJA is sup- 
ported by a Royal Society University Research Fellow- 
ship and by EPSRC under grants EP/I030298/1 and 
EP/EO30173. 
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SUPPLEMENTARY MATERIAL 



Proof of Eqs. (3) and (4) in the main text 



ated with the Malliavin weight. We therefore write 



We start by introducing the probability distribution 
P({r*i};i) for the particle positions, and the joint prob- 
ability distribution P({r*j}, q\; t) for the combination of 
the particle positions and the Malliavin weight. The two 
are related by 

P({?i};t) = fdq x P({r i },q x ;t). (SI) 

We also define the conjugate variable 

_ d]nP{{fi};t) 1 dP({n};t) 



and the conditional average 



jdqx q\P({n},q\;t) 
JdqxP({n},qx;t) 



(S2) 



(S3) 



The first goal is to show the equivalence of Q\ and 
(q\){ri} — Eq. (H]) in the main text. We do this by showing 
that they both obey the same evolution equation. 

As hinted at in the main text, we introduce an ex- 
plicit Euler-type scheme for updating the particle posi- 
tions. This also makes explicit the updating rule associ- 



DfiSt 



(S4) 



where r/ is the updated position of the ith particle at 
time t + St, St is the time step, and the & are a set of 
3iV Gaussian random variates of zero mean and variance 
2D5t. The corresponding updating rule for the Malliavin 
weight is 



, 1 J^dfo 

1\ = <?A + TTT^ > . — ' & 



2k B T f-f d\ 

2—1 



(S5) 



Note that the exact same sequence of random variates 
£j is used for updating the particle positions and the 
Malliavin weight. In this Euler scheme P({fi};t) obeys 
a Chapman-Kolmogorov equation 



P({r/}; t + 5t) = SUi d 3 rt P{{n}; t) 



x W({n'}\{n}) 



(S6) 



where the propagator is 
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W(m\m = (4ttJ> St)-™" exp (- <g ^ D. fl 5t/k B Tf j (g?) 

Differentiating the Chapman-Kolmogorov equation with respect to A leads to an adjunct equation for the conjugate 
variable Qx, 

Qx({n'};t + St)P({n'};t + 5t) = J ]Jd 3 n ^Q x ({f i }- t ) + ^^]p({f i };t)W({n'}\{n})^ (S8) 

i 

The second quantity in the square brackets in Eq. (|S8[) is 



9A <9A 2fc B P f-f dX V ' ' fc B T 



We now show that Eq. (|S8[) for the evolution of Q x is replicated by the evolution equation for the Malliavin weight. 
The joint probability distribution function P({r*i}, q x ; t) obeys an extended Chapman-Kolmogorov equation 

P({f/}, q' x ; t + 6t) = J [] d 3 n d 3 ^ dq x P({n}, q x ; t) P({£}) 5(r/ - n - A,) S(q' x - q x - A,) (S10) 

i 

where 

g) = ^ + g , A,(*,S) —Ef-Ci (SH) 

2 — 1 

give the discrete increments to the particle positions and Malliavin weight. The two 8- functions in Eq. (jSlOj) enforce 
these updating rules. The function P({£i}) is a 3iV-dimcnsional Gaussian. The 3iV random variates are uncorrelated 
and have zero mean and variance 2D5t. With the definitions in Eqs. (|S2j) and fjS3f) in hand we take the first moment 
of Eq. (|S10P with respect to q x to get 

(q\){? t >},t+6t P({n'h t + St) = fdq' x q' x P({n'}, q' x ; t + St) (S12a) 

= I Hi d 3 n d3 & d qx dq' x q' x P({ri},qy, t) P({£}) 5(r t > - n - A,) S(q' x -q x - A,) (S12b) 

= / Hi d 3 n d 3 £, dq x [q x + A q ] P({n}, q x ; t) P({£}) 5{r l ' - - A,) (S12c) 

= IHi d 3 n dq x [q x + A q ] P({r* }, qx] t) W{{n'}\{n}) (S12d) 

= m t d 3 n [(qx){? i} + WPdnhVWimm)- (S12e) 

To progress from Eq. (|S12b[) to (|S12e[) we do successively the q' x integral, the integrals, and the q x integral; using 
first the (^-functions then the definition of the conditional average for the final step. In doing the £j integrations, we 
recover the propagator W{r i '\ , ri) given by Eq. (|S7[) and the q x increment becomes 

i— 1 

This is identical to Eq. (fS9]) therefore we conclude that Eq. (|S12cj) for updating (qx)^ is identical to Eq. (fS8|) for 
updating the conjugate variable Q x . By choice these two quantities can be given the same initial values, thus 
establishing their equivalence, in other words 

JdqxqxP({fi},qx;t) 01nP({ft};f) 
jdq x P{{r l },q x -t) dX { * ' 

To establish Eq. in the main text we use Eq. (|S1[) and the second half of Eq. (|S2[) to rearrange Eq. (|S14|) into the 
alternative form / dq x q x P({r,}, g^; i) = <9P({r z }; t)/dX. Eq. © follows from this since 

(^9a) =JHi d 3 n dqx A({n}) qx P{{n }, g A ;t) 

rrr ^ Af^x^ ^MM dqn^ 3 ^({^})P({rl};0) _ d{A) (S15) 
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Practical implementation 

The implementation of MWS in a Brownian dynamics 
code is quite straightforward. The main point is that for 
efficiency one may want to update the Malliavin weight (s) 
according to Eq. (|S5|) at the same time as calculating the 
forces. This requires that the set of random variates £j 
be computed and stored at the start of the step, as in the 
following schematic algorithm : 

• generate the 37V random variates £j , 

• compute the forces fi{{fi}) and the quantity 
Etx(dfi/9X) ■ ii, 

• update the particle positions according in Eq. (|S4[) , 

• update the Malliavin weight (s) according in 
Eq. (J55]», 

• record any quantities of interest, i. e. and 
the value of the Malliavin weight (s). 

The last item does not necessarily have to be done every 
time step of course. As indicated in the main text, the 
algorithm is initialised by setting the particle positions to 
their initial values and the Malliavin weights to zero. For 
steady state problems there are two approaches to the use 
of the MWS correlation function method. The simplest 
is to choose a set of equally spaced reference points £o = 
nT for the calculation of Aq\ = q\{t) — qx(to), where 
T is a time period longer than the expected relaxation 



time of the system (which may have to be determined 
by trial and error). In this approach, every n timesteps 
the current values of q\{t) and A(t) are recorded. The 
running value of q\ is then set to zero and the simulation 
continued. The time average of the stored values of q\A 
gives d(A)/d\. More efficient is to use a sliding window 
to calculate the correlation function. Block averaging can 
be used for error estimates. 



Underdamped Brownian dynamics 

The extension of MWS to underdamped Brownian dy- 
namics is fairly simple. We denote the particle positions 
by ri and velocities by Uj. The dynamical equations are 

dri _ dvi - 

-jr = Vi , to— = fi - jVi + rji (S16) 

where to is the mass, 7 is the frictional drag coefficient, 
and ffi are white noise terms of amplitude 2^k^T . An 
Euler-type scheme for Eqs. (|S16[) is 

f/ = n+ ViSt Vi = Hi + — — h Ci (S17) 

to 

where the £j are 3N Gaussian random variates of zero 
mean and variance 2-/k B TSt/m 2 (note that we divided 
the velocity equations through by to). It follows that the 
probability distribution function P({rj}, {w;}; t) evolves 
according to the Chapman-Kolmogorov equation 



P({n'},{$iht + tt) = J n d3 ^ d3 ^ p ({^}'{^}; < ) 5 (^'-^-^^({^'}l{^}'{^}) ( S18 ) 

i 

where the partial propagator is 

W(m\Hln) = (^ferft/m»)-™/»^ . (S19) 

If we differentiate Eq. (|S18[) with respect to some parameter A we obtain 

Qx({n'}, {vi'}; t + 6t) P({n'}, {vi'}; t + St) 

= J Hi d 3 n d^lQ^in}, {Vi};t) + dlnW/dX] (S20) 

x s(n' - n - m) W({WWd, {fi}) P({fi}, {fi};t) 

1 

where Q\{{fi}, {fi}] t) = d In P({fi},{vi};t)/d\. This assuming the parameter of interest features only in the 

again suggests the updating rule for the Malliavin weight, force law, gives the Langevin equation 
q' x = q\ + dlnW/dX, in other words the derivation is 

impervious to the presence of a (5-function in the full , N „ 7 

propagator. Inserting the explicit expression for W, and — - = ■ — — — • ffi . (S21) 

dt 2jkBT ^ oX 

1=1 
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FIG. SI: (color) Numerical simulation of Eqs. (|S22|) and 
(fS23|) with k = 2, m = 1/2, x = 1, and 7 = T = 1. The 
falling curves show (x(t)) for ft = (solid line) and h — 0.1 
(dashed line). The rising curves show (x(t)qh(t)) (blue solid 
line) and d{x(t))/dh evaluated using forward finite differenc- 
ing (blue circles) . Averages are over 10 5 replicate simulations; 
cf. Fig. [1] in the main text. 



This result is very similar to the overdamped case. As 
a demonstration let us revisit the example of Brownian 
motion in a one-dimensional trap, but this time consider 
the underdamped case. The Langevin equations are 



dx 



dv 



— = v, ni- 
di dt 



and 



dqh 
dt 



27T 



(S22) 



(S23) 



Fig. IS II shows simulation results confirming that (qhx) 
d(x)/dh. 



Higher-order derivatives 

We next demonstrate how MWS extends to the compu- 
tation of higher-order derivatives. A double application 
of Eq. ([3]) in the main text, for two parameters A and /z, 
gives 



d 2 (A({n})) 
d\d[i 



(A({?i}) (qx„ + qxqj) 



(S24) 



where qx/i = dqx/dfi. Differentiating the discrete up- 
dating rule q' x = q\ + d In W j dX with respect to [i gives 



= + W/dXdfi. We insert the expression for 
W from Eq. (|S7[) into this, and simplify, to find the cor- 
responding Langevin equation 



dqx^ 
dt 



1 



N 

E 



2k&T f-fldXdu 11 k B T dX 

i — 1 



d_h 
dpi 



(S25) 



The new feature in this is a drift term (the last term) 
which has the consequence that (qxfj) = —{qxqfi}- This 
ensures that d 2 (A) / dXd n, = if A is a constant. 

Note that, as a result of the peculiar properties of 
the stochastic differential calculus, Eq. (|S25j) is not sim- 
ply found by differentiating Eq. ^ in the main text; 
rather one has to proceed via the discrete updating rules. 
Also note that the second-order Malliavin weight q\^ 
defined in Eq. (|S25[) must be combined with the two 
first-order Malliavin weights q\ and q^ to obtain the cor- 
rect weighted average in Eq. ([S24|) . We have tested the 
second-order MWS scheme for the trapped interacting 
particle cloud under shear, see for example Fig [3^ in the 
main text. 



One-dimensional trap 

Here we derive the expressions in Eqs. ([6]) for the tran- 
sient behaviour of a particle in a one-dimensional har- 
monic trap. Eqs. ([5]) in the main text are 



dx 

7/7 



2T 



These can be integrated to find 



it) =x e- Kt + -(l 

K 



) + / dt' e-^-t'^it') , 



q h (t) 



1 

2T 



dt'r)(t') 



(S27) 

Since x and qh are summed Gaussian random noises, it 
follows that P(x,qh',t) is a Gaussian — cf. §3.5.2 in The 
Theory of Polymer Dynamics M. Doi and S. F. Edwards 
(OUP, 1986). To characterise this Gaussian, it suffices 
to calculate the first and second moments. The first mo- 
ments follow immediately from Eqs. (|S27P and are (x) = 
xoe~ Kt + (h/K)(l — e~ Kt ) and (qh) = 0, as given in the first 
line of Eqs. ©. Therefore x - (x) = J*df e^*"*'^')- 
The second moment of x (the second line of Eqs. ([6])) is 



(x 2 ) - (x) 2 = / dt' J dt" e^C-*') e -«(*-*") x 2TS(t' - t") = |(1 - e~ 2Kt ) 



(S28) 
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where we have used (v(t')il(t")) = 2T5(t' — t"). Likewise the cross correlation term and the second moment of qh 
(the third line of Eqs. ((6])) are 

{*Qh) = t4 f dt' fdf'e-^-^ x 2T6(t' - t") = -(I - e~ Kt ) 

(S29) 



Two-dimensional trap in shear 

The quasi-FDT result in the main text follows from 
the steady state probability distribution P ss (x,y) for a 
particle in a two-dimensional trap under shear, which 
can be solved in closed form. Let us recall the Langevin 
equations for this problem, 



dx 



dy 



dt - KX + ™ + ^ dt 



-ny + T)y . 



(S30) 



where n is the trap strength and 7 is the shear rate. In 
common with the main text we set the particle mobility 
to unity and write temperature in terms of Boltzmann's 
constant, so that we can write D = T for the diffusion 
coefficient. 

From the Smoluchowski equation, the steady-state dis- 
tribution function corresponding to these Langevin equa- 
tions satisfies 



dJ x 
dx 



dJy 

dy 



= 



(S31) 



where the components of the probability current (flux) 
are 



J x = ( kx + 7?/)P ss - T 



dP ss 
dx ' 



Jy = -nyP ss - T 



dP ss 
dy 



(S32) 



Since the Langevin equations are linear, we expect that 
P ss will be a bivariate Gaussian so we write 



exp(- 



-\Ax 2 



.By 1 - Cxy) (S33) 



where A, B and C are coefficients, to be determined. The 
simplest way to proceed is to insert this as an ansatz into 
Eqs. (|S3ip and (|S32[) . to find that the coefficients have 
to satisfy 



2k = (A + B)T, 

An = {A 2 + C 2 )T, 

Bk = Cj+(B 2 + C 2 )T, 

Cn = Aj/2 + (A + B)CT. 



(S34) 



These four conditions arise from equating to zero the 
constant term and the coefficients of x 2 , y 2 and xy in 
Eq. (|S31j) . Although there are only three unknowns, 
Eqs. (|S34[) arc interdependent and admit the unique so- 
lution, 



A = 



4k 3 



C=- 



( 7 2 + 4k 2 )T' 
27K 2 



B = 



4/c 3 + 2k7 2 
(7 2 + 4k 2 )T ' 



(S35) 



( 7 2 + Ak 2 )T ' 

Hence the complete steady state distribution function is 



Pss = 



x cxp 



ttTV7 2 + 4k 2 

2k 3 x 2 + (2k 3 + K A j 2 )y 2 - 2jK 2 xy 



(S36) 



( 7 2 + 4k 2 )T 
For reference, the associated moments are 



T 

K 



(S37) 



Differentiating the steady state distribution with respect 
to the shear rate yields 



d\nP st 



$7 

2k 2 (2kx — 7J/) (72: + 2ny) 
( 7 2 + 4k 2 ) 2 T 



(S38) 



7 2 + 4k 2 



In the limit 7 — > this reduces to (qj)xy = xy/2T. An 
immediate application of this is to deduce the quasi-FDT 
result in the main text: 



d(xy) 



d~i 



7=0 



<*V> 

2T 



(S39) 



